/*
 * EBFV1_calculateNonvisc.cpp
 *
 *  Created on: 14/09/2012
 *      Author: rogsoares
 */

#include "EBFV1_hyperbolic.h"

namespace PRS{

	double EBFV1_hyperbolic::calculateIntegralAdvectiveTerm(int dom, int dim, TMesh* pTMesh_data){
		double FluxIJ[3], alpha_max;

		//ofstream fid;fid.open("nonvisc.txt");

		int nedges,n;
		pTMesh_data->edge_list[dom].getSize(nedges,n);
		for (int i=0; i<nedges; i++){
			EdgeData* edata = pTMesh_data->edge_list[dom].getValue(i,0);
			VertexData* vdata_0 = edata->vertex_0;
			VertexData* vdata_1 = edata->vertex_1;
			double Sw_I = vdata_0->Sw;
			double Sw_J = vdata_1->Sw;

			// get high order approximation for saturation
			if ( pSimPar->useHOApproximation() ){
				pHOApproximation->getSw_HighOrderApproximation(pSimPar,edata,dim,Sw_I,Sw_J);
			}

			double fwII = pPPData->getFractionalFlux(Sw_I);
			double fwJJ = pPPData->getFractionalFlux(Sw_J);
			double fwIJ = 0.5*(fwII + fwJJ);

			// Numerical Flux Function
			for(int j=0;j<dim;j++){
				FluxIJ[j] = fwIJ*edata->vel_new[j];
			}

			// Fractional Flux Flow Function Derivative (with respect to saturation)
			double df_dsIJ = ( fabs(Sw_I-Sw_J) > 1.e-12 )?fabs((fwJJ-fwII)/(Sw_J-Sw_I)) : .0;

			// Approximate Eigenvalue (Note that we are using the linearized form of df_dsIJ)
			double alpha = norm_L2(edata->vel_new,dim)*df_dsIJ;

			// get the maximum alpha to compute the time step
			//if ( alpha > alpha_max) alpha_max = alpha;
			alpha_max = std::max(alpha,alpha_max);

			// Central difference Contribution
			double non_visc_fv = dot(FluxIJ,edata->Cij,dim);

			// Numerical Diffusion
			double non_visc_ad = 0.5*edata->Cij_norm*alpha*(Sw_J - Sw_I);

			// Computing "Non-Viscous" Terms
			vdata_0->nonvisc = vdata_0->nonvisc + (non_visc_fv - non_visc_ad);
			vdata_1->nonvisc = vdata_1->nonvisc - (non_visc_fv - non_visc_ad);
		}
	//	fid.close();
	//	STOP();


		// todo: find a MPI function alpha_max = P_getMaxDbl(alpha_max);
		double dt_new = (pSimPar->CFL()*pTMesh_data->getMaximumEdgeLentgh()*pSimPar->getPorosity(3300))/alpha_max;
		double dt_old = this->getTimeStep();
		this->setTimeStep( std::min(dt_new,dt_old) );
	}

	void EBFV1_hyperbolic::resetNodalNonviscTerms(int dom, TMesh* pTMesh_data){
		int nrows, ncols;
		pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
		for (int j=0; j<nrows; j++){
			pTMesh_data->vertices_list[dom].getSize(nrows,ncols);
			VertexData* vdata = pTMesh_data->vertices_list[dom].getValue(j,0);
			vdata->nonvisc = .0;
		}
	}
}
